#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <mpi.h>

double objs[] = {
  0. ,0., -100.5, 10000., 0., 0., 0., 0.25,
  0.272166, 0.272166, 0.544331, .027777,
  0.643951, 0.172546, 0., .027777,
  0.172546, 0.643951, 0., .027777,
 -0.371785, 0.099620, 0.544331, .027777,
 -0.471405, 0.471405, 0., .027777,
 -0.643951,-0.172546, 0., .027777,
  0.099620,-0.371785, 0.544331, .027777,
 -0.172546,-0.643951, 0., .027777,
  0.471405,-0.471405, 0., .027777,
  4., 3., 2., 1., -4., 4., -3., 1., 5.
};

double *intersect(double x, double y, double z, double dx, double dy, double dz, double *maxp)
{
  int i;
  double *o = objs, *oo = 0;
  double max = *maxp;
  double xx, yy, zz, b, t;

  for (i = 0; i < 11; i++)
    {
      xx = *o++ - x; 
      yy = *o++ - y; 
      zz = *o++ - z;
      b = xx * dx + yy * dy + zz * dz;
      if ((t = b * b - xx * xx - yy * yy - zz * zz + *o++) < 0 || (t = b-sqrt(t)) < 1e-6 || t > max)
	continue;
      oo = o - 4;
      max = t;
    }
  *maxp = max;
  return oo;
}

double shade(double x, double y, double z, double dx, double dy, double dz, int de)
{
  double max = 1e6, c = 0, r, k, *o;
  double nx, ny, nz, ldx, ldy, ldz, rdx, rdy, rdz,inv;
  int i;

  if (!(o = intersect(x, y, z, dx, dy, dz, &max)))
    return 0;
  x += max * dx; y += max * dy; z += max * dz;
  nx = x - *o++; ny = y - *o++; nz = z - *o++;
  r = sqrt(nx * nx + ny * ny + nz * nz);
  //inverse of r only calculated once
  inv = 1.0/r;
  nx *= inv; 
  ny *= inv; 
  nz *= inv; 
  k = nx * dx + ny * dy + nz * dz;
  rdx = dx - 2 * k * nx; rdy = dy - 2 * k * ny; rdz = dz - 2 * k * nz;
  o = objs + 44;
  for (i = 0; i < 3; i++)
    {
      ldx = *o++ - x; 
      ldy = *o++ - y; 
      ldz = *o++ - z;
      r = sqrt(ldx * ldx + ldy * ldy + ldz * ldz);
      inv = 1.0/r;
      ldx *= inv; 
      ldy *= inv; 
      ldz *= inv; 
      if (intersect(x, y, z, ldx, ldy, ldz, &r))
	continue;
      if ((r = ldx * nx + ldy * ny + ldz * nz) < 0)
	continue;
      c += r;
      //power function replaced
      if ((r = rdx * ldx + rdy * ldy + rdz * ldz) > 0)
      {
       double r_2 = r*r*r;
       double r_3 = r_2 * r_2;
       r_3 *= r_3;
       r_3 *= r_2;
       c += 2 * r_3;
      }
	 }
  if (de < 10)
    c += .5 * shade(x, y, z, rdx, rdy, rdz, de + 1);
  return c;
}


void calc_tile(int size, int xstart, int ystart, int tilesize, char* tile)
{
  double dx, dy, dz, c, r;
  int x, y, i;
  //inverse outside loop
  double inv = 1.0/(float)(size-1);
  i=0;

  for (y = ystart; y < ystart+tilesize; y++)
    for (x = xstart; x < xstart+tilesize; x++)
      {
        
        double xx = x *inv;
        double yy = 1. - y * inv;
        dx = -0.847569 - xx * 1.30741 - yy * 1.19745;
        dy = -1.98535  + xx * 2.11197 - yy * 0.741279;
        dz = -2.72303                 + yy * 2.04606;
        //changed to inverse of r
        r = 1.0/sqrt(dx * dx + dy * dy + dz * dz);
	c = 100 * shade(2.1, 1.3, 1.7, dx * r, dy * r, dz * r, 0);
	if (c < 0)
	  c = 0;
	if (c > 255)
	  c = 255;
	tile[i++]=c;
      }
}

int main(int argc, char** argv)
{
  int size=4800,tilesize=100,xtiles,ytiles;
  double total_time=0;
  /* number of tiles in x and y direction */
  xtiles=size/tilesize;
  ytiles=size/tilesize;

  MPI_Init(NULL, NULL);  
  int rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  int world_size;
  MPI_Comm_size(MPI_COMM_WORLD, &world_size);
  
  //Each of them define the tile
  char *tile;
  if((tile=(char*)malloc(tilesize*tilesize*sizeof(char)))==NULL)
  {
     fprintf(stderr,"Could not allocate tile memory!\n");
     MPI_Abort(MPI_COMM_WORLD, 1);
  }

  int xleft=rank*xtiles/world_size;
  int xright=xleft+xtiles/world_size;
  int i,xc,yc;
  
  if(rank==0)
  {
    char *picture;
    int count=0;
    int tilebase;
    //Define the picture only in the master process
    if((picture=(char*)malloc(size*size*sizeof(char)))==NULL)
    {
      fprintf(stderr,"Could not allocate picture memory!\n");
      MPI_Abort(MPI_COMM_WORLD, 1);
    }
    for(i=0; i<size*size*sizeof(char); ++i)
     picture[i]=0;
    
    //modify directly first part (only in the rank 0)
    double t1 = MPI_Wtime();
    for(yc=xleft; yc<xright; yc++)
    {
      for(xc=0; xc<ytiles; xc++)
      {
	/* calc one tile */
	calc_tile(size, xc*tilesize, yc*tilesize, tilesize, tile);
	tilebase=yc*tilesize*tilesize*xtiles+xc*tilesize;
	for(i=0; i<tilesize; i++)
	{
	    memcpy((void*)(picture+tilebase+i*tilesize*xtiles),
		   (void*)(tile+i*tilesize),
		   tilesize*sizeof(char));
	}
	count++;
      }
    }
    
    double t2=(MPI_Wtime()-t1);
    total_time+=t2;

    int j;
    t1 = MPI_Wtime();
    for(j=1;j<world_size;++j)
    {
      //int myCon=0;
     for(yc=j*xtiles/world_size; yc<j*xtiles/world_size+xtiles/world_size; yc++)
     {
       //printf("Tile %d\n",yc);
       for(xc=0; xc<ytiles; xc++)
	 {
	   /* Recieve the tile (it is not calculated)*/
	   tilebase=yc*tilesize*tilesize*xtiles+xc*tilesize;
	   MPI_Recv(tile,tilesize*tilesize, MPI_CHAR , j , 0 , MPI_COMM_WORLD , MPI_STATUS_IGNORE);
	   //printf("Recieved from %d \n",j);
	   for(i=0; i<tilesize; i++)
	   {
	    memcpy((void*)(picture+tilebase+i*tilesize*xtiles),
		   (void*)(tile+i*tilesize),
		   tilesize*sizeof(char));
	   }
	   count++;
	   //myCon++;
	   //printf("\nCounter receptor %d\n",myCon);
	 }
     }
     //     printf("\nProcess: %d\n",j);
    }
    t2=(MPI_Wtime()-t1);
    total_time+=t2;
    
    // UNCOMMENT FOR WRITING BITMAP

  FILE *fd=fopen("result.pnm","w");
  fprintf(fd,"P5\n%d %d\n255\n",size,size);
  for(count=0; count<size*size; count++)
    fputc(picture[count],fd);
  fclose(fd);

    free(tile);
    free(picture);
  }
  else
  {
    printf("\nProcess: %d\n",rank);
    int con=0;
     double t1 = MPI_Wtime();
     for(yc=xleft; yc<xright; yc++)
     {
       for(xc=0; xc<ytiles; xc++)
	 {
	   /* calc one tile and send it*/
	   calc_tile(size, xc*tilesize, yc*tilesize, tilesize, tile);
	   //printf("Sent from %d \n",rank);
	   MPI_Send(tile,tilesize*tilesize,MPI_CHAR,0,0,MPI_COMM_WORLD);
	   con++;
	   //printf("\nCounter sender %d\n",con);
	 }
     }     
     double t2=(MPI_Wtime()-t1);
     total_time+=t2;
     free(tile);
  }

  MPI_Finalize();
  printf("\n Time: %f\n",total_time);
  exit(0);
}
